PROBABILISTIC REPRESENTATION OF BERNOULLI, EULER, AND 
CARLITZ HERMITE POLYNOMIALS 



C. VIGNAT, UNIVERSITE D'ORSAY 



Abstract. We revisit in a probabilistic framework the umbral approach of Bernoulli numbers, 
Euler numbers and Carlitz Hermite polynomials by Gessel pp. This study allows to explicit 
equivalents of some famous umbras. 



1. Introduction 

In PP, Gessel shows how the umbral calculus, introduced by Blissard and later developed by Rota 
and Roman [2J, allows to derive some elementary but also more elaborate results about Bernoulli 
numbers and about classical polynomials such as Charlier and Hermite polynomials. In this note, 
we show that another approach based on probabilistic representations can be substituted to the 
umbral formalism. This substitution brings interesting results about the relationship between 
random variables and umbrae. 



2. Bernoulli numbers and Bernoulli polynomials 
The sequence of Bernoulli numbers {B n } is defined by its generating function 



t n t 



5^ s »-r = -r-r' 1*1 < 277 

^— ' n\ e t — 1 

n=0 

and the sequence of Bernoulli polynomials {B n (x)} by their generating function 

+°° te^ x 

V B n (x) — = -, \t\ <2tt : x£ R, 

* — ' n\ e l — 1 

so that the Bernoulli numbers verify B n — B n (0) . First values are 

1 1 3 x 

B (s) = 1; B x (x)=x- -; B 2 (x) = x 2 - x + -; B 3 (x) = x 3 - -x 2 + - 

and 

11 1 

Bq = 1: Bi = ; Bo = — ; B-\ = 0; Ba = . 

In [3], Sun gives the following probabilistic representation of the Bernoulli polynomials 

Theorem. [Sun] Given a sequence {L n } of independent random variables, each with Laplace 
distribution i exp (— |x|) , x € M, define the random variable 

+oo T 

(2.1) L = yp~. 

fe=i 

Then the following probabilistic representations hold^ 

(2.2) B n {x) = E (iL + x --\ , n > 0, x £ R 



2 



and 



(•'_>.:{) B„ = E ( iL- i ) . /, > n. 



^In this paper, the symbol E denotes the expectation operator Eg (X) = f f (X) g (X) dX where / is the 
probability density of the relevant random variable. 
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The random variable L, being denned as a series of i.i.d. random variables, is not easy to 
characterize. Thus we propose the equivalent result. 

Theorem. The random variable in (|2.1I) follows a logistic distribution, with density 

fh (x) — — sech (nx) , iGR. 
Proof. The random variable L in (|2.ip has characteristic function 

E (e ltL ) = L 

v 1 sinh(|) 

But from 1.9.2] 

f + °° l 2 / \ / \ j TTt , (lXt 

I sech (ax) cos (xt) ax — — r-cscn — 
Jo 2 « 2 V 2a 

so that, with a = tt, the density of L is 



/l (x) = -sech (ttx) . 



We note from [5J p. 471] that the random variable L can be obtained as 



□ 



L = — log = — log — 

2tt 8 1 - U 2ir & E 2 

where U is uniformly distributed on [—1, +1], E\ and E 2 are independent with exponential distri- 
bution /e (x) — exp (—x) , x G [0, +oo[ and equality is in the sense of distributions. 

Corollary 1. The Bernoulli polynomials read 

(2.4) B n (x) = E(— log — — + x - =e(— log — +x-- 

K ' y> \2tt h l-U 2 J \2n & E 2 2 

and the Bernoulli numbers 

B n = Ef^log^- 1 -) 71 = E (^\o g ^- 1 -Y , n>0. 

\2tt h l-U 2 J \2tt h E 2 2) ~ 

We remark that the result by Sun was already given by Talacko in [3] . 
We note that L is also [9] a Gaussian scale mixture 

L = 2K.N 

where N is Gaussian and K follows the Kolmogorov-Smirnov distribution 

fx (a?) = 8x^2 n 2 exp (-2n 2 x 2 ) , x > 0. 

n=l 

This result was given first by Barndorff- Nielsen et al. in [TU], where the mixing distribution of 2K 
is also characterized by its moment generating function 



ip K (s) = E 



sin (tv\/2s) 

Proposition 2. The Bernoulli numbers satisfy 

B n =E(iL+^\ , n^l. 
Proof. We first compute the generating function 

+oo n +oo , 
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as it should be. Now let us consider 

+00 / 

El 1 \ i iC I 



l\ n t n te l t 



n=0 



so that, since E(zL + 3) = 5, 



x ^ s n ^ n +00 , \\ n t r 



which shows the result. □ 

Let us show now that the representation (|2.3[) allows to recover easily some fundamental results 
about Bernoulli numbers. 

In umbral calculus, as described in [lj, a Bernoulli number B n is represented by an umbra B 
which should be replaced by B n each time it appears as B n . 

Proposition 3. [Gessel, (7.2)] The Bernoulli numbers satisfy, in the umbral notation, 

(B + l) n = B n , Vn^l. 

Proof. In a probabilistic setting, this is a direct consequence of the result of Proposition [5] for 
n + 1, 

B n = E ( iL - \ ) = E (iL + i 



Proposition 4. The Bernoulli numbers verify, in the umbral notation, 

(B + l) n = {-l) n B n , n^l. 

Proof. Since by symmetry 

f 1 E x 1\" ( 1 E 2 1 

B — log-i-- = E — log-f-- 



□ 



2?r £: 2 2 / V 2vr Si 2 



(-ir*(s">«|+5)"-<- 1 »"* 



for n > 2, we deduce the result. □ 

Other integral representations of Bernoulli numbers and polynomials exist in the litterature; for 
example, (|2.4H can be considered as a substitute for the Mellin Barnes integral [T51 24.7.11] 

1 r-c+100 2 

B n (x) = — / (x + t) n ( —\ dt; 

moreover, the integral representation [6, p. 39, eq. (27)] 

r+00 

B 2n = (-1)" +1 7T / t 2 "csch 2 (nt) dt 
Jo 



can be easily deduced from (|2.3I) ; however, this identity holds only for even order Bernoulli numbers 
- the odd order ones being equal to except B\ , as a consequence of propositions [3] and 01 This 
can be a difficulty in the calculation of some series which involve all Bernoulli numbers. As an 
example, we provide a short proof of Kaneko's theorem. 

Theorem. [Kanekoj The Bernoulli numbers verify 

I n . ) (n + i + 1) B n+l = 0, n > 0. 

i=0 ^ ' 
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Proof. It can be easily checked that Vz € C, n > 0, 

n+l 



J2[ i + i + z n +'' = (n + 1) z" (1 + z) n (1 + 2z) 



so that 



n+l 

E 

i=0 



n+l 



(n + i + l)B n+l = (n + l)E[iL 



iL 



(2iL) 



2i (n + l)E 



L —L 



The random variable L being symmetric, i.e. —L being distributed as +L, the former expectation 
can be computed as 



E 



L\-L z - 



1 



= E 



-L -L 



1 



= 



what proves the result 

Kaneko's theorem is a special case of Momiyama's identity 
'm + 1 

j (n + H + l)U n+k = {-L)- - 

k=0 



(-i) m E 



k=0 



k 



B n+k - Y.( n l l ) (m + k+ 1) B m+k 



which can also be easily proved by remarking that 



E 

fe=0 



m + 1 
k 



{n + k + l)B n+k = E L {n+l)[iL-~) [iL+- 



1 



(to + 1) [iL - - [%L + - 



2 

n+l 



m+1 



1 



(to + n + 2) ( iL - - 



m+n+l 



which, by the symmetry of L, can be shown to coincide with the right-hand side sum. 
Another famous identity can be proved easily using the probabilistic approach 



Theorem 5. The following identity holds 

B n+i (y) - B n+ i (x) 



rv 

/ B n (z) dz = 

J X 



n + l 



, n > 2 



Proof. Using the probabilistic representation (|2.2[) . 



J X 



z) dz 



E 



1 



n + 1 



y / i 

lL + z ~2 



E[iL + y 



n+l 



E { iL + x 



B n+1 (y) - B n+1 (x) 
n + l 



An equivalent version of Theorem [S] is 

-^B n (z) = nB n -i (z) 

which shows that Bernoulli polynomials are Appell polynomials. 
A related result is the following [13 24.13.2] 

Theorem 6. The Bernoulli polynomials satisfy 

nX+l 



(2.5) 



B n (z) dz = x r ' 



□ 



□ 
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Proof. From Thm [5l this integral is equal to 

B n+1 (x + 1) - B n+1 (g) _ E (iL + x + \) n+1 - E (iL + x - \) n+1 
n + 1 n + 1 

Expanding each power with the binomial formula yields 



— Y 

+ 1 ^ 

k=0 



n+1 / , -, \ ( / -, \ n+1 — k / \ n+1— fc ' 



By Q, all terms vanish except for k = n, which yields the result. □ 

Another identity, used by Gessel in the proof of Kaneko's formula, and also proved in |7| using 
the extended Zeilberger's algorithm, reads 

Theorem 7. The Bernoulli numbers satisfy the following identity 

^2 f )B n+k = (-l) m+ "^ ( k )B m +k. 

k=0 ^ ' fc=0 ^ ' 



Proof. The proof is straightforward using (|2.2I) since the left-hand side is 



while the right-hand side is E (iL — i) (iL + ■§) • Since L is distributed as —L, 

which concludes the proof. □ 
In [7], K.W. Chen considers the sequence of numbers {K n } , n > defined as 

= ^ ( . )-B n+ i_|_i, 



i=0 v 7 

and proves the following 

Theorem 8. [Chen] The sequence {K n } satisfies the identities 

<A (2n - k\ 2n 

M * / 2n~k = ~ 2 "' n - 1 

and 



2-4 /2n - fc - l\ 2n-l 

fc=0 v 7 



Proof. We give a short proof of the first identity remarking that K n — E i(iL — \) n+1 (iL + |) ' 
and that 



so that 



fc=0 
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with U = iL - \ and V = %L + \ so that UV = —L 2 — | and y/l - AUV = iL. Finally, 



E 

fc=0 



2n - k\ 2n 



k J 2n — k 



K k = E iL- 



-iL 



2n 



iL 



EUL-- 



2n+l 



1 



2 



1 



The first term is -E?2n+1 — while the second reads 

2 n 



E[iL--\E[iL 



1 



E\iL--\E\iL+~ 



2n 



E iL 



1 \E \iL 



i i + 1 



B, 



k-i- 



Bn ~ o(l-o») ^' 
A quick proof can be given using the following 
Lemma 10. TTie power sum polynomial has integral representation v7c,Vn £ N, 

S k (n) = / S fc (z + 1) jfc. 



1 n— 1 / \ 
_£«*( W_ fc («-1 



./o 

Proof. This can be verified considering n £ K in the right-hand side of (|2.6|) so that 

k / 1 \ / -, \ k — i 



l^n) - a£(*)»<<-l>"(.L-i 



= # L I + n + - 



= E L [iL + n 



and, since SV- (0) = 0, 



S k (n) 



and the result holds in particular Vn £ 



E L [iL + z 



dz = B n (z + 1) dz 



Another simple proof of this lemma can obviously be deduced from identity (|2.5p . 
Using this integral representation, we compute now 



n — l / \ n / \ 

Y,^[l)BkS n - k {a-l) = ^ k [ n k )BkS, 

k=0 ^ ' k=0 ^ ' 



(a -I)- a n (a - 1) B„ 



The right-hand side sum is 



k=0 



iL\ — 



1 



k r a-l 



B n - k (z + l)dz= E Lu l 



a iL\ — 



□ 



which is equal to i?2n+i — Bin — —Bin since n > 1. The proof of the second identity is equally 
easy. 

A last and less easy identity is the following from 
Theorem 9. [Gesselj Denote the power sum polynomial 

(2.6) ft(») = E < * = E(- 1 ) fc ^ /M " 

t=l i=0 

Then, with a £ N, a > 1 and n£N, 



□ 



«L 2 + z 



which can be integrated as 



1 



n + 1 



Bl u l 2 



1 



a iLi — — — 1L2 + a — — 



1 



n+l 



a I iL\ — - I — I iL 2 + i 



1 



n+l 
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Expanding both (n + 1) powers, we obtain 

g (" 1 >"«-' ((* + - i)' - - + f " 

But by ([3]), all terms cancel in this sum except for k = 1 and k — n in which cases they add up, 
so that we obtain 

— !— (- (n + 1) a n B n ) + — j— ((n + 1) ofl„) = B„a (l - a"- 1 ) 
n + 1 n + 1 

which, adding the remaining term —a™ (a — 1) B n , yields the result. 

3. EULER NUMBERS AND EULER POLYNOMIALS 

3.1. definition and characterisation. We derive here analogous results for Euler numbers and 
Euler polynomials. A generating function for the sequence {E n } of Euler numbers is 

^ En ni =sech w 

ra=0 



and for the Euler poynomials {E n (x)} 

+ 

The Euler numbers are obtained as 



— t n 2e tx 



n\ e l + 1 

n— 



E n = 7T~ E n ( — 

2™ \ 2 



First values are 



1 3 1 

E (x) = 1; Ex (x)=x- -; E 2 (x) = x 2 - x; E 3 (x) = x 3 - -x 2 + - 

and 

E = 1; E x = 0; E 2 = -1; E 3 = 0; E A = 5. 
In [3], the following formula is derived 

Theorem 11. [Sun] If Lq is defined as 

+ oo , 

where are independent and Laplace distributed, then the Euler polynomials read 
(3.2) E n {x) = E (iL + x-^ 

and the Euler numbers 

E n =2 n E n (^j = 2 n E(iL ) n . 

We provide a more convenient characterization of the random variable Lq as follows. 
Theorem 12. The random variable Lq follows the hyperbolic secant distribution 

f Lo (x) = sech (ttx) . 
Proof. The characteristic function of Lq is 

Ee lLot = sech ( l 



From [3 1.9.1], 

f +aD 7T / 7T 

/ sech (ax) cos (xt) dx — — sech ( — ; 
J y ' v ' 2a \2a 

so that, with a — ir, the density of Lq is 

f Lo (x) = sech (ttx) . 
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Thus ttLq follows an hyperbolic secant distribution. □ 

We note from [5] that the random variable Lq can be obtained as 
(3.3) L = -log|C| = i(Iog|JVi|- log|JV 2 |) 

7T 7T 

where C is Cauchy distributed and N\ and N2 are two independent standard Gaussian random 
variables. 

The random variable Lq is also a scale mixture of Gaussian with mixing distribution given in 
|10) . the moment generating function of which reads 

<P ( s ) 



cos (tt\/2s) 



At last, the random variable iLq has the same moments as the Levy area - that is the signed area 
- of a Brownian motion B t for < t < 1 , see [H] . 

We leave to the reader the proofs of basic identities such as 

(3.4) E n (l-x) = (-l) n E n (x), 

( U )E r {x) = E n (x + l) 

r=0 ^ 7 

and 

(3.5) E( 2 2 "V = ° 

r=0 ^ 7 

as a consequence of the probabilistic representation (|3.2[) : the first identity is a consequence of the 
symmetry of the distribution of Lq, the second is obtained using the binomial theorem. 

3.2. links between Bernoulli and Euler polynomials. 

3.2.1. a summation identity. An interesting identity that links Bernoulli and Euler polynomials is 
the following [13, p. xxxiii] 



1 " / \ 
B n {x) ^ — Y J [i)Bn-kE k {2x 



2 

k=0 



and its more general version [T5l 24.14.5] 



k=0 



The proof reads 

1 A/n\ If J 1 \ " f L + L e x + y 1 

J B ^W^W = ^ liL e + x--+iL + y--) = E ( 1 ^— H — — — — 



2 n ^ \kj w v 7 2 n V 2 2/ V 2 2 2 

fe=0 v 7 v 7 v 

and to conclude, we need simply to show that Lo ~^ Le ~ L e : this can be deduced immediately from 
the identities (|3.1[) and (|2.1[) ; or from 

£o + £ e 1 , , 1 , fW 1, IATxIVET 



2 2?r 1 1 2tt V E 2 2?r ° V^l^; 

JV 

distributed, we deduce 



= — log|C| + — log- 



where C is Cauchy, E is exponential and since -J= is Laplace distributed and \L\ is exponentially 



1 1 £l r 
— log — ~ L P . 

2n & E 2 
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3.2.2. an integral identity. Another similar identity is |15l 24.13.3] 
(3.6) ( X+ ~ 2 B n {z)dz- En[2x) 



2 n+l 



which can be proved as follows 



J ' B n (z) dz = — I E (iL e + x) n+1 - E (^iL e + x 



1 \^(,L e L Q 2x x " +1 / T T 



E i^r- + i-r- + -TT- 





- 1 ! 








1 


'.Le 


H 


La 






2 


2x - 


2 


)) 



n+llV2 2 27 V 2 2 2 

2» +1 FTT) { E {{ iL ° + 1) + { iLo + 2x - CO"" - E - 1) + + 2x " { 

By property all terms cancel in the binomial expansions of the (n + 1) powers, except for k = n 
and the only remaining term is 

1 fn+l\ ( 1\" E n (2x) 

! )E [iLo + 2x - - = 



2»+! (n + 1) V n / V V 2 ™ +1 ' 

3.2.3. a consequence of (|3.6p . As a consequence of the identity (|3.6[) . we deduce the identity 

E n (x + l) + E n (x) = 2x n 

as follows: 

E n (x + l) + E n (x) = 2 n+1 I B n (z) dz = 2 n+1 (-J = 2a;' 1 
where we applied the result (|2.5[) . 

4. HERMITE POLYNOMIALS 

In [T], the study of Hermite polynomials involves the definition of the umbra M such that 

e Mx = e -x 2 ^ xe£ ^ 

In the rest of this section, we provide a probabilistic interpretation to this umbra, and using a 
master identity on Gaussian random variables, derive new simple proofs of the results by Gessel. 

Theorem 13. The umbra M is equivalent to the expectation 

(4.1) f(Mx) = Ef(-iZx) 

where Z is a Gaussian random variable with zero mean and variance a 2 — 2, which we denote as 
Z ~ A/"(0, 2), for any admissible function f. 

Proof. It suffices to prove 

Ee~ iZx = e-* 2 

which is nothing but the expression of the characteristic function of the Gaussian random variable 
Z. □ 

We need the following master identity. 

Proposition 14. In umbral notation, 

(A 0) Mx+M 2 y _ 1 „ XD ( ^_ 

Proof. Expressing 

e Mx+M 2 y _ E Mx+MVy _ fiyg-fr+Vy) 

where V ~ AT (0, 2) , we need to compute the expectation 

1 



e -^ e -i-+yyf dV . 

oc 
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We recognize here, up to a constant term ^L, the convolution, evaluated at the point -, of a 
variable V ~ Af (0, 2) with a random variable W ~ Af ^0, ^p-^j ; it can easily be checked to be the 
Gaussian with variance N (0, 1 + 4j/ 2 ) , i.e. 



: exp 



y/i + Ay 2 V 1 + 4 2/ 



As special cases, we recover 

(1) for y = 0, e Mx = exp (-a; 2 ) 

(2) for x = 0, = ^ 

Here are now a few results from pQ that can be deduced from the master identity (|4.2D . 

(1) the umbral notation for the Hermite polynomials 

H n (u) = (2u + M) n 
can be deduced from the identity (|4.1D with / (M) = (2u + M) n , so that 

H n (u) =E{2u + iN) n 

where N ~ Af (0, 2) , which is a well-known result (see for example [331 P-49]) 

(2) the generating function 

1 / 4m 2 x 



Ei?2n (u) — r = ; CX P i 

n— v 



can be derived from 



tec n / 2\ 

H 2n (u) = E N exp (Ax (u + iN) 2 ^ = E v exp ( - (v^fx - iuV^x) ) 

n=0 ' ^ 7 

and application of the master identity (|4.2[) 
(3) the same approach yields the bivariate generating function 



Ex y 6 4x 

w m! n! VI + 4x 



//•),, T "_ 

2y^ 



2 > 



V 



1 + Ax 



5. Carlitz and Zeilberger's Hermite polynomials 
In pQ , Gessel proposes an umbral study of Carlitz Hermite polynomials 



in An 



lu ..i k y n k 



,.k 



□ 



fc=0 v 7 v 7 

and of the Zeilberger's Hermite polynomials 

mAn , \ / \ 

^w=g(r)(;)^ 

with notation m An = min (m, n) . In this aim, he defines the umbrae A and B by 

(5.1) A m B n = 5 m , n ml, Vm, neN 
or equivalently as 

(5.2) exp (Ax + By) = exp (xy) Vs, y e R. 

The probabilistic counterpart of these umbrae is given by the following result. 
Proposition 15. Let us consider a complex circular normal random variable Z with density 

/ z (z) = iexp(-|z| 2 ). 

The umbras A and B in (|5.1I) and (|5.2p are identified with the expectations with respect to Z and 
Z respectively. 
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Proof. It suffices to check that 

EZ m Z n = 6 m , n m\ 

which is straightforward. 
Equivalently, 

E exp (Zx + Zy) = E exp {(X + iY) x + (X — iY) y) 

= E exp (X (x + y)) E exp (iY (x - y)) 

= exp^(x + y) 2 cxp-^(x-y) 2 

so that the result holds. □ 

As a consequence, we deduce a probabilistic representation of Zeilberger's Hermite polynomials 

as 

H„ hn (u)=E{l + Z) m (l + uZ) n 

and of the Carlitz Hermite polynomials as 

H m , n (u, v) = u m ~ n E (1 + Z) m (uv + Z) n . 

Using now classical results about complex Gaussian random variables, we can recover easily some 
of the results derived by Gessel using umbral calculus. We need the following master identity, the 
umbral version of which is Lemma 5.3 in pQ. 

Lemma 16. If Z is circular normal, then 
(5.3) 

Eexp (Za + Zb + Z 2 u + ZZv + Z 2 w) = 1 exp ( ^ + b 2 u + ab (I - v) \ 

^J(l-v) 2 -4uw \ t 1 -^ ~ Auw J 

Proof. The proof is a consequence of the fact that the expectation in (|5.3p can be rewritten as 

E exp [Za + Zbj 

where Z is a complex Gaussian variable with covariance matrix 

1 



EZZ* = E 



\J (1 — v) 2 — 4xu 



x 

(l-v) 2 



The latest expectation over Z is nothing but the generating function of Z computed at the point 
(a, b). □ 

In order to illustrate the efficiency of the probabilistic approach, we provide some quick deriva- 
tions of the results obtained by Gessel. 

(1) the generating function of Zeilberger's Hermite polynomials reads 

V H min (w)^-r^-r = E z exp(x(l + Z))exp(y(l + wZ)) 
L — ' m! n\ 

m,n-0 

= exp (x + y) Ez (xZ + ywZ) 
= exp (x + y + wxy) 

where the latest equality is a consequnce of (|5.2[) 

(2) the bilinear generating function of Zeilberger's Hermite polynomials reads 

I OO tl 

V H mtn (u)H mjn (v)^-^- = Eexpx(l + Z 1 )(l + Z 2 )exp(y(l + uZ 1 )(l + vZ 2 )) 
z — / ml n\ 

m,n— 

= E exp (x (1 + Z 2 ) + y (l + Z 2 v)) 

x exp (xZi (1 + Z 2 ) + uyZ x (l + vZ 2 )) 

Taking first expectation on Z\ and using formula (|5.2I) . then expectation on Z 2 and using 
formula (|5.3p yields the result. 
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(3) the following general formula holds 



(5.4) 




i\ jl kl 11 ml ~~ ^ 



v l U! J x k y l t 




(1 — ut) 2 — 4:U 2 xy 



1 



x exp 



(1 + uw) 2 x + (1 + uv) 2 y + Auxy + (1 — ut) (v + w + t + uvw) 



(1 — ut) 2 — 4u 2 xy 



The quintuple sum can be computed as 



E'exp v (1 + Z) + w (l + uZ) + x (1 + Z) 2 + y (l + uZ) 2 +t (1 + Z) (l + uZ) 



and using the master lemma (|5.3I) and simple algebra, we recover (|5.4p . 



6. Conclusion 



We have seen that all umbrse used in [T] have a simple probabilistic expectation counterpart. 
Depending on the complexity of the formula to prove, each approach has its advantages and 
drawbacks. A future direction of work is to deepen the understanding of the link between the 
umbral and the probabilistic approach. 



[1] I.M. Gessel, Applications of the classical umbral calculus, Algebra univers. 49 (2003) 397-434 
[2] S. Roman, The Umbral Calculus, Academic Press, 1984 

[3] P. Sun, Moment representation of Bernoulli polynomial, Euler polynomial and Gegenbauer polynomials, Sta- 
tistics & Probability Letters 77 (2007) 748-751 

[4] J. Talacko, A note about a Family of Perks' Distribution, Sankhya, The Indian Journal of Statistics, 20-3-4, 
323-328, Dec. 1958 

[5] H. Bateman, Tables of Integral Transforms, vol.1, McGraw-Hill Book Company, 1954 
[6] H. Bateman, Higher Transcendental Functions, vol.1, Krieger Pub Co, 1981 

[7] K.-W. Chen, A summation on Bernoulli numbers, Journal of Number Theory 111 (2005) 372-391 
[8] L. Devroye, Non-Uniform Random Variate Generation, Springer, 1986 

[9] L. A. Stefanski, A normal scale mixture representation of the logistic distribution, Statistics & Probability 
Letters 11 (1990) 69-70 January 1991 

[10] O. Barndorff-Nielsen, J. Kent and M. Sorensen, Normal Variance-Mean Mixtures and z Distributions, Interna- 
tional Statistical Review, 50-2, 145-159, Aug. 1982 

[11] I. Gessel, Solution to problem E3237, Amer. Math. Monthly 96 (1989) 364 

[12] D. Levin and M. Wildon, A combinatorial method for calculating the moments of Levy area, Trans. Amer. 

Math. Soc. 360 (2008) 6695-6709 
[13] I. S. Gradshteyn, I.M. Ryzhij, Table of integrals, series, and products, Academic Press, Seventh Edition, 2007 
[14] S. Janson, Gaussian Hilbert Spaces, Cambridge University Press, 1997 

[15] F.W.J. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clarck, NIST Handbook of Mathematical Functions, NIST 
and Cambridge University Press, 2010 



References 



